% !TeX root =thesis.tex
\chapter{The variational approach  using  the BCS-type ansatz\label{ch:mean}}
We can also formulate the problem  using the  BCS-type ansatz.  It is not difficult to calculate  the expectation of the free energy as well as  the wave function that optimizes the free energy within this ansatz.  The optimization process gives us gap equations, which determine the  wave function with the constraint from number equations.   %with those common parameters as chemical potential, gap.  
We will see that this method yields  the mean-field solution consistent with the one of the path-integral approach.  However, this method is mean-field by nature and  difficult  to be extended for studying collective excitations.  

We rewrite the Hamiltonian Eq. \ref{eq:pathInt2:ham2} in momentum space, and restore the momentum dependence of the interaction coefficients.  
\begin{equation}\label{eq:uvw:hamiltonian}
\begin{split}
 H=&\sum_\vk\epsilon^a_\vk{}a^+_\vk{}a^{}_\vk+\sum_\vk\epsilon^b_\vk{}b^+_\vk{}b^{}_\vk+\sum_\vk\epsilon^c_\vk{}c^+_\vk{}c^{}_\vk\\
  &-\sum_{\vk\vk'}U_{\vk\vk'}a^+_\vk{}b^+_{-\vk}{}b^{}_{-\vk'}a^{}_{\vk'}
	-\sum_{\vk\vk'}V_{\vk\vk'}a^+_\vk{}c^+_{-\vk}{}c^{}_{-\vk'}a^{}_{\vk'}\\
 &-\sum_{\vk\vk'}Y_{\vk\vk'}a^+_\vk{}b^+_{-\vk}{}c^{}_{-\vk'}a^{}_{\vk'}
	-\sum_{\vk\vk'}Y^*_{\vk\vk'}a^+_{\vk'}{}c^+_{-\vk'}{}b^{}_{-\vk}a^{}_{\vk}
\end{split} 
\end{equation}
Here we only keep the zero center-of-mass momentum paring terms as in the original BCS work\cite{BCS}. 
We take the free atom at zero magnetic field as the  zero energy.  
\begin{equation*}
\epsilon_{\vk}^i=k^2/(2m)+\eta^i, \qquad (i=a,b,c)
\end{equation*}
 $\eta^i$ is the Zeeman energy of the $i^{\text{th}}$ atom at magnetic field $B$. The hermiticity  of the Hamiltonian imposes
 \begin{equation}
 U_{\vk'\vk}=U^*_{\vk\vk'},\qquad{} V_{\vk'\vk}=V^*_{\vk\vk'}
\end{equation}
  We then introduce the BCS-type ansatz 
\begin{equation}\label{eq:ansatz}
 \ket{\Psi}=\prod_\vk\br{u_\vk+v_\vk{}a^\dg_\vk{}b^\dg_{-\vk}+w_\vk{}a^\dg_\vk{}c^\dg_{-\vk}}\ket{0}
\end{equation}
$\ket{0}$ is the particle vacuum state.  We require $\abs{u_\vk}^2+\abs{v_\vk}^2+\abs{w_\vk}^2=1$ for normalization.  This ansatz is constructed like the original BCS ansatz for superconductivity.  An alternative ansatz would be $\prod_\vk(u_\vk+v_\vk{}a^\dg_\vk{}b^\dg_{-\vk})(u^{'}_\vk+w_\vk{}a^\dg_\vk{}c^\dg_{-\vk})\ket{0}$, which is actually the same as Eq. \ref{eq:ansatz}, because  the cross term vanishes due to the Pauli exclusion on the common species.   As the original BCS ansatz, this ansatz does not have the fixed particle number.  We  will just require the expected value of number operators match the total particle number. For all  interaction terms, we get  two types of contributions,
for example, 
\begin{equation*}
\av{U_{\vk\vk'}a^\dg_\vk{}b^\dg_{-\vk}{}b^{}_{-\vk'}a^{}_{\vk'}}
=\sum_{\vk}U_{\vk\vk}\abs{v_\vk}^2+\sum_{\vk\neq\vk'}U_{\vk\vk'}v^{}_{\vk'}u^*_{\vk'}u^{}_\vk{}v^*_\vk
\end{equation*}
The first term is the Hatree term and the second term is the more interesting pairing term.  


This gives the full free energy as 
\begin{equation}\label{eq:uvw:F}
 \begin{split}
  &F\equiv\av{H-\mu{}N}\\
    =&\sum(\xi^a_\vk+\xi^b_\vk)\abs{v_\vk}^2+\sum(\xi^a_\vk+\xi^c_\vk)\abs{w_\vk}^2\\
    &-\sum_{\vk}U_{\vk\vk}\abs{v_\vk}^2-\sum_{\vk\neq\vk'}U_{\vk\vk'}v^{}_{\vk'}u^*_{\vk'}u^{}_\vk{}v^*_\vk\\
    &-\sum_{\vk}V_{\vk\vk}\abs{w_\vk}^2
      -\sum_{\vk\neq\vk'}V_{\vk\vk'}w^{}_{\vk'}u^*_{\vk'}u^{}_\vk{}w^*_\vk\\
    &-\sum_{\vk}Y_{\vk\vk}w^{}_{\vk}v^*_\vk{}
      -\sum_{\vk\neq\vk'}Y_{\vk\vk'}w^{}_{\vk'}{u^{*}_{\vk'}}v^*_\vk{}u^{}_\vk\\
    &-\sum_{\vk}Y^*_{\vk\vk}w^*_{\vk}v^{}_{\vk}{}
      -\sum_{\vk\neq\vk'}Y^*_{\vk\vk'}w^*_{\vk}{u^{}_{\vk}}v^{}_{\vk'}{}u^{*}_{\vk'}
 \end{split}
\end{equation}
where we have set
\begin{equation*}
 \xi^a_\vk=\epsilon^a_\vk-\mu^a,\qquad\xi^b_\vk=\epsilon^b_\vk-\mu^b,\qquad\xi^c_\vk=\epsilon^c_\vk-\mu^b
\end{equation*}
 Two chemical potentials are added to make sure the $n_a=n_b+n_c=\nth{2}n$.  In principle, $\mu^{a}$ does not need to be equal to $\mu^{b}$  as there is no exchange or conversion between $(a)$ and $(b,c)$.  In fact, the structure of the ansatz guarantees that $n_a=n_b+n_c$. Therefore, we set $\mu^{a}=\mu^{b}$ for simplicity and drop the superscript on chemical potential hereinafter. 
We will also drop the Hartree terms because they are  only related to the density and thus can be absorbed into the chemical potentials.   In the second summation, we will ignore the fact that the summation only goes through $\vk\neq\vk'$ because the corrections due to this restriction lead to only the higher order terms. 
 
 \section{Exact gap equations and number equations}
We introduce two new parameters $h_{1\vk}$ and $h_{2\vk}$, which corresponds to the mean field values (equal time average) of the anomalous Green's functions,  (in the same way as  $h_{1,2}$  in Sec \ref{sec:pathInt2:meanfield}), 
\begin{gather}
u_{\vk}^2+v^{2}_{\vk}+w^{2}_{\vk}=1\\
u_{\vk}v_{\vk}=h_{1\vk}\\
u_{\vk}w_{\vk}=h_{2\vk}
\end{gather}
We can solve $u_{\vk}$, $v_{\vk}$, $w_{\vk}$ in terms of  $h_{1\vk}$ and $h_{2\vk}$ (all parameters are taken as real)\footnote{It is not hard to prove that the parameters that optimize the free energy are real, within an overall phase factor, when the interactions are real. (See also footnote \ref{foot:pathInt2:real} in Chapter \ref{ch:path2}, page \pageref{foot:pathInt2:real}.)}.  One complication in solving above equations is that $u_\vk$, $v_\vk$ and  $w_\vk$ are all  monotonic functions of $\vk$, while $h_{1\vk}$ or $h_{2\vk}$ is not in the BCS end.  So, we need to be careful when taking the square root.  We introduce $\sgn_k$ for such purpose.  
\begin{equation}\label{eq:uvw:uvwh12}
\begin{split}
u_{\vk}^2=&\frac{1}{2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)\\
v^{2}_{\vk}=&\frac{2 h_{1\vk}^2}{1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}
=\frac{ h_{1\vk}^2}{2( h_{1\vk}^2+ h_{2\vk}^2)} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)\\\
w^{2}_{\vk}=&\frac{2 h_{2\vk}^2}{1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}
=\frac{h_{2\vk}^2}{2( h_{1\vk}^2+ h_{2\vk}^2)} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)
\end{split}
\end{equation}
  where $\sgn_k=1$  for all $k$ in BEC cases, and  for large $k$ in BCS cases. In BCS cases,  $\sgn_{k}=-1$ when $k$ is small.  In the single-channel crossover, $\sgn_{k}=\sgn(\epsilon_{k}-\mu)$, and the turning point is at $u_\vk^2=1/2$ which corresponds to zero chemical potential, $\mu=0$.  The two-channel crossover is more delicate to treat.  The turning point is still at $u_\vk^2=1/2$, but it is no longer exactly at $\mu=0$.  The $\sgn_{k}$ function  is very important in number equations Eqs. \ref{eq:20100909:number}.  They however can be  absorbed later on, after we convert the formulas back into notations of  $(u_{\vk},\,v_{\vk}, \,w_{\vk})$ or $\xi_{k}$.


%Below, we adopt $\sgn_k=1$. (This works for BEC and most part of BCS; and in the final equations, the sign function is taken care of automatically by using $E_{\vk}$.  It is easy to verify that  we arrive at the same final result in the other cases where $\sgn_{k}=-1$. ) 
Eq. \ref{eq:uvw:uvwh12}'s  derivatives over $h_{1\vk}$ are\footnote{Here we follow the same convention by taking $h_{1\vk}$ and $h_{2\vk}$ as real}
\begin{equation}
\begin{split}
\pdiff{u_{\vk}^2}{h_{1\vk}}=&-\frac{2 h_{1\vk}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\sgn_{k}\\
\pdiff{v_{\vk}^2}{h_{1\vk}}=&\frac{2 h_{1\vk}\sgn_k}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}-\frac{8 h_{1\vk} h_{2\vk}^2\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\\
\pdiff{w_{\vk}^2}{h_{1\vk}}=&\frac{8 h_{1\vk} h_{2\vk}^2\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\\
=&\sgn_{k}\frac{h_{1\vk} h_{2\vk}^2 \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}{2\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} (h_{1\vk}^2 + h_{2\vk}^2)^2}
\end{split}
\end{equation}
The second equation is not very obvious, but can be obtained by noting that $\pdiff{v_{\vk}^2}{h_{1\vk}}=-\pdiff{u_{\vk}^2}{h_{1\vk}}-\pdiff{w_{\vk}^2}{h_{1\vk}}$.   Similarly, derivatives over $h_{2\vk}$ are
\begin{equation}
\begin{split}
\pdiff{u_{\vk}^2}{h_{2\vk}}=&-\frac{2 h_{2\vk}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\sgn_{k}\\
\pdiff{v_{\vk}^2}{h_{2\vk}}=&\frac{8 h_{1\vk} ^{2}h_{2\vk}\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\\
=&\sgn_{k}\frac{h_{1\vk}^{2}h_{2\vk} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}{2\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} (h_{1\vk}^2 + h_{2\vk}^2)^2}\\
\pdiff{w_{\vk}^2}{h_{2\vk}}=&\frac{2 h_{2\vk}\sgn_k}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}-\frac{8 h_{1\vk} ^{2}h_{2\vk}\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}
\end{split}
\end{equation}
We can obtain the gap equations by differentiating the free energy with respect to $h_{1\vk}$ and $h_{2\vk}$.
\begin{subequations}\label{eq:20100909:fullgap}
\begin{multline}
\frac{h_{1\vk}\sgn_k}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}} \xi^{ab}_{\vk}+\frac{4 h_{1\vk} h_{2\vk}^2\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\eta\\
-\sum_{\vk'}U_{\vk\vk'}h_{1\vk'}-\sum_{\vk'}Y_{\vk\vk'}h_{2\vk'}=0
\label{eq:20100909:fullgapa}
\end{multline}
\begin{multline}
\frac{h_{2\vk}\sgn_k}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}} \xi^{ac}_{\vk}-\frac{4 h_{1\vk}^{2} h_{2\vk}\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2} \left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\eta\\
-\sum_{\vk'}V_{\vk\vk'}h_{2\vk'}-\sum_{\vk'}Y_{\vk\vk'}h_{1\vk'}=0
\label{eq:20100909:fullgapb}
\end{multline}
\end{subequations}
where $\eta=\epsilon^{ac}_{\vk}-\epsilon^{ab}_{\vk}=\eta^{c}-\eta^{b}$ is the bare Zeeman energy difference and is larger than most other energy scales, such as $E_{F}$.  It is close to  the binding energy of the closed-channel bound state when it is not too far away from the resonance point.   We see one  disadvantages of the ansatz method: unlike the path integral approach, where $h_{1\vk}$ and $h_{2\vk}$ have a very specific forms, Eqs. (\ref{eq:pathInt2:h1}, \ref{eq:pathInt2:h2}), and ready for further approximation and analysis,  we do not have much guide for these two quantities here and need to make some guess.  

We can define order parameters
\begin{gather}
\Delta_{1\vk}=\sum_{\vk'}U_{\vk\vk'}h_{1\vk'}+\sum_{\vk'}Y_{\vk\vk'}h_{2\vk'}\\
\Delta_{2\vk}=\sum_{\vk'}V_{\vk\vk'}h_{2\vk'}+\sum_{\vk'}Y_{\vk\vk'}h_{1\vk'}
\end{gather}
Both of them should have little $\vk$-dependence at low momentum for short-range interactions.  We will drop the $\vk$ subscripts in both $\Delta$ in the following. 
Two number equations can be obtained.  One for each channel.
\begin{subequations}\label{eq:20100909:number}
\begin{gather}
N_{\text{open}}=2\sum_{\vk}v_{\vk}^{2}
	=\sum_{\vk}\frac{ h_{1\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)
	\label{eq:meanfield:opennum}\\
N_{\text{close}}=2\sum_{\vk}w_{\vk}^{2}
	=\sum_{\vk}\frac{ h_{2\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)
\end{gather} 
\end{subequations}


\section{Approximate solution of the mean field equations}
Here instead of solving these equations, we just demonstrate that they are consistent to the saddle point solution obtained in the path-integral method up to the first order of correction $\zeta$.

Let us look at the closed-channel gap equation (Eq. \ref{eq:20100909:fullgapb}) first.  It can be rewritten as 
\begin{multline*}
\frac{h_{2\vk} \sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\br{2 \xi_{\vk}+\eta-\frac{4 h_{1\vk}^{2}}{\left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\eta}
-\sum_{\vk'}V_{\vk\vk'}h_{2\vk'}-\sum_{\vk'}Y_{\vk\vk'}h_{1\vk'}=0
\end{multline*}
We can put some terms back into the notation of $(u_{\vk},\,v_{\vk},\, w_{\vk})$.
\begin{gather}
\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}=(u^{2}_{\vk}-v_{\vk}^{2}-w_{\vk}^{2})\\
1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}=2u_{\vk}^{2}\label{eq:meanfield:1sqrt}\\
\frac{4 h_{1\vk}^{2} }{\left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}=\frac{v_{\vk}^{2}}{u_{\vk}^{2}}
\end{gather}
We rewrite the closed-channel gap equation using these relations
\begin{equation*}
\frac{h_{2\vk}}{u^{2}_{\vk}-v_{\vk}^{2}-w_{\vk}^{2}}\br{2 \xi_{\vk}+\eta-\frac{v_{\vk}^{2}}{u_{\vk}^{2}}\eta}-\sum_{\vk'}V_{\vk\vk'}h_{2\vk'}-\sum_{\vk'}Y_{\vk\vk'}h_{1\vk'}=0
\end{equation*}
At high momentum ($\gg{}k_{F}$),  the third term in the parenthesis is negligible and the the denominator $u^{2}_{\vk}-v_{\vk}^{2}-w_{\vk}^{2}\approx1$.  This equation is then very similar to the two-body \sch equation (Eq. \ref{eq:intro:close}). For the closed-channel, the interaction, or the closed-channel bound state, $\phi_{0}$, has non-zero value for a range  in the momentum space much larger than $k_{F}$.  Thus we expect the closed-channel correlation $h_{2\vk}\propto\phi_{0}f(\vk)$ where $f(\vk)$ approach one at high momentum.  This is in fact the same argument as Eq. \ref{eq:pathInt2:hphif}  in the beginning of Chapter \ref{ch:path2}.  The next thing to notice is that the third term in the parenthesis might dominate the other two terms in  certain situation (low momentum, BCS) because the denominator  $u_{\vk}^{2}$ can be very small in such regions.  A nature way to compensate this big term is to set\footnote{$u_{\vk}^{2}$ in Chapter \ref{ch:path2} is defined to be $(\xi_{\vk}+E_{\vk})/2E_{\vk}$ and does not carry an obvious physical meaning as $u_{\vk}^{2}$ in this appendix. They are the same in the lowest order of our expansion over $\zeta$ for the narrow resonance. This does not affect our analysis because the closed-channel correlation is actually a first order quantity.} $f(\vk)=u_{\vk}^{2}$.  Now we arrive the same conclusion for the closed-channel correlation as in the path-integral approach (Eq. \ref{eq:pathInt2:hphi})
\begin{equation}\label{eq:meanfield:h2phi}
h_{2\vk}=\tilde\alpha{}\phi_{0\vk}u_{\vk}^{2}
\end{equation}
Note that $\tilde\alpha$ is not fixed yet and should not be identified as $\alpha$ in the path-integral approach (Eq. \ref{eq:pathInt2:hphi}) \textit{a prior}. The equation becomes,
\begin{equation*}
\frac{\tilde\alpha{}\phi_{0\vk}u_{\vk}^{2} }{u^{2}_{\vk}-v_{\vk}^{2}-w_{\vk}^{2}}\br{2 \xi_{\vk}+\eta-\frac{v_{\vk}^{2}}{u_{\vk}^{2}}\eta}-\sum_{\vk'}V_{\vk\vk'}\tilde\alpha{}\phi_{0\vk}u_{\vk}^{2}-\sum_{\vk'}Y_{\vk\vk'}h_{1\vk'}=0
\end{equation*}
Using the similar approach as in the path-integral method, we multiply the equation with $\phi_{0\vk}^{*}$ and integrate the equation in order to find $\tilde{\alpha}$. In the denominator of the first term $u^{2}_{\vk}\gg{}v_{\vk}^{2}$, $u^{2}_{\vk}\gg{}w_{\vk}^{2}$ in the most of the integral domains (up to order of $\kappa$($\eta$) in the momentum space).  The equation is approximately 
\begin{equation*}
\tilde\alpha{}\sum_{\vk}\frac{\phi_{0\vk}\phi_{0\vk}^{*}u_{\vk}^{2}}{u^{2}_{\vk}}\br{2 \xi_{\vk}+\eta}-\tilde\alpha\sum_{\vk\vk'}V_{\vk\vk'}{}\phi_{0\vk}\phi_{0\vk}^{*}u_{\vk}^{2}-\sum_{\vk\vk'}\phi_{0\vk}^{*}Y_{\vk\vk'}h_{1\vk'}=0
\end{equation*}
And it is not hard to find 
\begin{equation}
\tilde\alpha\approx\frac{\sum_{\vk\vk'}{\phi_{\vk}^{*}}{Y_{\vk\vk'}}{h_{1\vk'}}}{{-E_{b}+\eta-2\mu-\tilde\lambda_{1}}}
\end{equation}
 And $\tilde\lambda_{1}$ is simpler than Eq. \ref{eq:pathInt2:lambda1}
\begin{equation*}
\tilde\lambda_{1}(\eta)\equiv-\sum_{\vk\vk'}{\phi_{\vk}^{*}}{v_{\vk'}^{2}V_{\vk\vk'}}\phi_{\vk'}
\end{equation*}
But it should be equal to $\lambda_{1}$ from the path-integral approach at the first order of $\zeta$ (Eq. \ref{eq:pathInt2:lambda1}) because this term is the dominated one (see Sec. \ref{sec:pathInt2:lambda}, Eq. \ref{eq:pathInt2:lambda1es}). We will just use $\lambda_{1}$ hereinafter. $\tilde\alpha$ is the same as $\alpha$ in the path-integral approach if $h_{1\vk}$ here is the same as the one in the path-integral approach,  which we will show in the following. Particularly for low momentum, $\phi_{0\vk}\sim\nth{\eta+\epsilon_{\vk}}$ (see Appendix \ref{sec:pathInt2:short-range}) and it is not hard to find at low momentum
\begin{equation}\label{eq:meanfield:aphi}
\alpha\phi_{0\vk}\overset{{k\lesssim{k_F}}}{=}\frac{\Delta_2}{\eta}
\end{equation}

   Now let us check the open-channel gap equation.  We can rewrite it as 
   \begin{equation*}
   \frac{h_{1\vk} \sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\br{ \xi^{ab}_{\vk}+\frac{4 h_{2\vk}^2}{\left(1+\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2}\eta}=\Delta_{1}
   \end{equation*}
Again, using Eq. \ref{eq:meanfield:1sqrt}, we replace the denominator with $2u_\vk^2$.  Also we replace $h_{2\vk}$ with $\tilde\alpha{}\phi_{0\vk}u_{\vk}^{2}$ using Eq. \ref{eq:meanfield:h2phi}.  
 \begin{equation*}
   \frac{h_{1\vk}\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\br{ 2\xi_{\vk}+({\tilde\alpha}\phi_{0\vk})^2\eta}=\Delta_{1}
   \end{equation*}
 At low momentum, the second term in the parenthesis can be simplified as $\Delta_2^2/\eta=\zeta\Delta_{1}$  using Eq. \ref{eq:meanfield:aphi}
   \begin{equation*}
   \frac{h_{1\vk}\sgn_{k}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}\br{ 2\xi_{\vk}+\zeta\Delta_{1}}=\Delta_{1}
   \end{equation*}
   We can solve $h_{1\vk}$ in terms of $\Delta_1$.  
   \begin{equation*}
   h_{1\vk}^2=\frac{\Delta_1^2(1-4h_{2\vk}^2)}{(2\xi_{\vk}+\zeta\Delta_{1})^2+4\Delta_1^2}
   \approx\frac{\Delta_1^2(1-4h_{2\vk}^2)}{4E_\vk^2+4\xi_\vk\Delta_{1}\zeta}
   \approx\frac{\Delta_1^2}{4E_\vk^2}(1-\frac{\xi_\vk\Delta_{1}}{E_\vk^2}\zeta)(1-4\frac{\Delta_{1}}{\eta}u_\vk^4\zeta)
   \end{equation*}
   Here we take the low momentum value of $h_{2\vk}^{2}\approx\frac{\Delta_{1}}{\eta}u_\vk^4\zeta$ in the last step. This term is much smaller  than $\frac{\xi_\vk\zeta}{E_\vk^2}$ because of the factor $\frac{\Delta_{1}}{\eta}$ at the low momentum.  We can neglect this term. Now we have 
    \begin{equation}\label{eq:meanfield:h1approx}
   h_{1\vk}^2\approx\frac{\Delta_1^2}{4E_\vk^2}(1-\frac{\xi_\vk\Delta_{1}}{E_\vk^2}\zeta)
   \end{equation}
   This is equivalent in the first order of $\zeta$ to the $h_{1\vk}$ in the path-integral approach (Eq. \ref{eq:pathInt2:h1}). 
   
   
It is easy to see that $h_{2\vk}^2$ is the higher order as $\zeta$ and $w_\vk\ll1$ all the time.  So we always have 
$h_{2\vk}\approx{}w_\vk$ except in low momentum.  However, the summation in the low momentum is very small, one order higher in $\zeta$.  The closed-channel number equation is then simply $N_{\text{closed}}=\sum{}w_\vk^2\approx{}\sum_{\vk}h_{2\vk}^2$, which is the same result in the path-integral approach as we discussed in Eq. \ref{eq:pathInt2:closeD2}.

We move to the open-channel number equation Eq. \ref{eq:meanfield:opennum}.  At the momentum smaller than the characteristic momentum of the closed-channel bound state $\kappa$  ($\eta\sim{}E_{b}=\kappa^{2}/2m$), factor $\frac{ h_{1\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)} $ is approximately $1$,
\begin{equation*}
\frac{ h_{1\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)} =1-\frac{h_{2\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)}
\overset{k\ll{}\kappa}\approx1-\frac{h_{2\vk}^2}{h_{1\vk}^2}
\approx1-\frac{4E_{\vk}^{2}}{\eta\Delta_{1}}u_{\vk^{4}}\zeta
\sim1+O(\zeta^{2})
\end{equation*}
while on the high momentum 
\begin{equation*}
\frac{ h_{1\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)}
\overset{k\gtrsim\kappa}\approx\frac{\frac{\Delta_{1}^{2}}{4E_{\vk}^{2}}}{\frac{\Delta_{2}^{2}}{\eta+\epsilon_{\vk}}}
\sim\frac{\Delta_{1}^{2}}{\epsilon_{\vk}\Delta_{2}}\ll1
\end{equation*}
Therefore the factor $\frac{ h_{1\vk}^2}{( h_{1\vk}^2+ h_{2\vk}^2)} $ limit the summation into  low momentum.  On the other hand, in the low momentum we have, $\frac{h_{2\vk}^{2}}{h_{1\vk}^{2}}\sim\frac{\Delta_{1}}{\eta}\zeta$, thus we can neglect $h_{2\vk}^{2}$ in the square root comparing to $h_{1\vk}^{2}$. And using Eq. \ref{eq:meanfield:h1approx}
\begin{equation*}
\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\approx\sqrt{1-4h_{1\vk}^{2}+O({\zeta}^{2})}
	\approx\frac{\xi_{\vk}}{E_{\vk}}+\frac{\Delta_{1}^{2}}{4E_{\vk}^{3}}\zeta+O({\zeta}^{2})
\end{equation*}
So we have the open-channel number equation 
\begin{equation}
N_{\text{open}}\approx\sum_{\vk}1-\frac{\xi_{\vk}}{E_{\vk}}-\frac{\Delta_{1}^{3}}{4E_{\vk}^{3}}\zeta
\end{equation}
This is consistent to the open-channel number equation derived from the path-integral approach  (Eq. \ref{eq:pathInt2:openNum}, where the second term is negligible comparing to the third term).

Thus we have shown that the mean-field solution of the path-integral approach can solve the gap equations and number equations of the ansatz method up to the first order of $\zeta$. 
%At the lowest order, the second term in both gap equations, Eq. \ref{eq:20100909:fullgap}, are negligible.  And to the lowest order, the first term in these two equations are equal $\Delta_{1}$ and $\Delta_{2}$ respectively.  











































%The number equation is (see footnote (\ref{foot:20100909:sgn}) for $\sgn_{k}$)
%\begin{equation}\label{eq:20100909:number}
%N=2\sum_{\vk}(v_{\vk}^{2}+w_{\vk}^{2})=\sum_{\vk}\left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)
%\end{equation} 
%where N is the number of atoms (twice  the number of pairs)
%
%\section{Approximating closed-channel with bound-state level}
%Let us look at \eef{eq:20100909:fullgapb}.  It is energetically very disadvantageous to deviate from the bound state as the binding energy and the difference from the next bound state is much  larger than Fermi energy.     This equation must close to  the solution for the isolated closed-channel \sch equation.  
%\begin{equation}\label{eq:20100915:twobody}
%{\phi^{0}_{\vk}}(2\epsilon^{ac}_{\vk})-V_{\vk\vk'}\phi^{0}_{\vk'}=-E_{b}\phi^{0}_{\vk'}
%\end{equation}
%This equation relates region (in k-space) much larger than $k_{F}$.  On the other hand, the open-channel 
%$h_{1\vk}$ is only substantial in region around or not too high above  $k_{F}$.  Therefore, as the first order approximation, we takes $h_{1\vk}=0$ for this equation. So we can drop the second term and the denominator in the first term ($h_{2\vk}\ll1$ all the time) of \eef{eq:20100909:fullgapb}\footnote{ This is certainly OK for he BEC end where $F_{k}\ll1$ all the time.  It is less satisfactory, might be problematic,  when $F_{k}$ is close to maximum $\nth2$, this happens around Fermi energy in BCS limit, but the approximation is still OK for other places. So at least for bulk of the region of $h_{2\vk}$ satisfies \eef{eq:20100915:gapb}}.  In  open-channel equation, \eef{eq:20100909:fullgapa}, we drop the factor $\left(1+\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)^2$ in the second term for easier calculation, (no strong reason yet, but this term can only takes value between 1 and 4, and the second term is relatively minor comparing to the first). We write down the approximated gap equations
%\begin{subequations}\label{eq:20100915:gap}
%\begin{gather}\label{eq:20100915:gapa}
%\frac{2h_{1\vk}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}} (\epsilon^{ab}_{\vk}-2\mu+  h_{2\vk}^2\eta)-U_{\vk\vk'}h_{1\vk'}-Y_{\vk\vk'}h_{2\vk'}=0\\
%\label{eq:20100915:gapb}
%{2h_{2\vk}}(\epsilon^{ab}_{\vk}-2\mu+\eta)-V_{\vk\vk'}h_{2\vk'}-Y_{\vk\vk'}h_{1\vk'}=0
%\end{gather}
%\end{subequations}
%Here we express $\xi^{ab}_{\vk}=\epsilon^{ab}_{\vk}-2\mu$, $\epsilon^{ac}_{\vk}=\epsilon^{ab}_{\vk}+\eta$.  And many-body closed-channel wave function $h_{2\vk}$ is proportional to two-body bound-state wave function $\phi_{\vk}^{0}$,  $h_{2\vk}=\alpha\phi^{0}_{\vk}$.  Comparing to $h_{2\vk}$ in path integral approach (Eq. \ref{eq:pathInt2:hphi} in Sec. \ref{sec:pathInt2:ren}), the many-body effect factor $u_{\vk}^{2}$ is missing.  It can be derived if taken into account the factor $\nth{{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}}}$ with more care.  However, contrast with a known structure of $h_{2\vk}$ there, we are quite blank with the behavior of $h_{2\vk}$, only in hinder sight, we realize this extra many-body factor, $u_{\vk}^{2}$.   Here we just stick to the lowest order and ignore this factor.  
%\eef{eq:20100915:gapb} becomes
%\begin{equation}\label{eq:20100915:GF}
%h_{2\vk}=-\frac{Y_{\vk\vk'}h_{1\vk'}}{2(E^{0}-\eta+2\mu)}
%\end{equation}
%Plug this back into the last term of \eef{eq:20100915:gapa}, we have 
%\begin{equation}
%\frac{2h_{1\vk}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}} (\xi^{ab}_{\vk}+  h_{2\vk}^2\eta)-U_{\vk\vk'}h_{1\vk'}+\frac{Y_{\vk\vk'}Y_{\vk'\vk''}}{2(E^{0}-\eta+2\mu)}h_{1\vk''}=0
%\end{equation}
%Now considering the last two terms has weak dependency in low k, we can set 
%\begin{equation}\label{eq:20100915:gap1}
%\frac{2h_{1\vk}}{\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}} (\xi^{ab}_{\vk}+  h_{2\vk}^2\eta)\equiv\Delta_{\vk}=(U_{\vk\vk'}h_{1\vk'}-\frac{Y_{\vk\vk'}Y_{\vk'\vk''}}{2(E^{0}-\eta+2\mu)}h_{1\vk''})
%\end{equation}
%We can express $h_{1\vk}$ according to $\Delta_{\vk}$,  (ignore the higher order of $h_{2\vk}$)
%\begin{equation}\label{eq:20100915:F}
%h_{1\vk}=\frac{\Delta_{\vk}}2\sqrt{\frac{(1-4h_{2\vk}^{2})}{(\xi^{ab}_{\vk}+  h_{2\vk}^2\eta)^{2}+\Delta_{\vk}^{2}}}
%\end{equation}
%Put it back into the second half of gap equation (\ref{eq:20100915:gap1}), 
%\begin{equation}\label{eq:20100915:onechannel}
%\Delta_{\vk}=\sum_{\vk'}\br{U_{\vk\vk'}-\frac{Y_{\vk\vk''}Y_{\vk''\vk'}}{2(E^{0}-\eta+2\mu)}}\frac{\Delta_{\vk'}}2\sqrt{\frac{(1-4h_{2\vk'}^{2})}{(\xi^{ab}_{\vk'}+  h_{2\vk'}^2\eta)^{2}+\Delta_{\vk'}^{2}}}
%\end{equation}
%At low k, $\Delta_{\vk}$ has weak dependency on k and we can take $\Delta_{\vk'}$ out of the summation,  then go through the same procedure as in \cite{Leggett,Fetter} to renormalize the equations. 
%
%\section{Number equation}
%For number equation \eef{eq:20100909:number}, more approximation can be made in Feshbach resonance.  closed-channel component $h_{2\vk}$ is much extended in k-space due to the tight-binding, and $h_{1\vk}$ is significant up to the order of Fermi energy $E_{F}$.  It seems OK to assume that most weight of closed-channel lies beyond $E_{F}$, so when $\epsilon_{k}\gg{E_{F}}$, the integrand becomes simply $h_{2\vk}^{2}$, the summation gives $N_{close}$.  Therefore we have 
%\begin{equation}\label{eq:20101004:number}
%N_{open}=\sum_{\vk}{}^{'}\left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)-h_{2\vk}^{2}
%\end{equation} 
%Here the summation only goes up to the order of $E_{F}$.
%
%
%Put \eef{eq:20100915:F} into number equation \eef{eq:20100909:number}, we have\footnote{Here we drop $\sgn_{k}$ as $(\epsilon^{ab}_{\vk}-2\mu+  h_{2,0}^2\eta)$ takes care of the sign}
%\begin{equation}\label{eq:20101011:number}
%N=\sum{1-\sqrt{1-4h_{2\vk}^{2}}\frac{(\epsilon^{ab}_{\vk}-2\mu+  h_{2\vk}^2\eta)}{{\sqrt{{(\epsilon^{ab}_{\vk}-2\mu+  h_{2\vk}^2\eta)^{2}+\Delta^{2}}}}}}
%\end{equation}
%This equation also requires to take into account that $h_{2\vk}$ approaches 0 at high momentum to be convergent.
%\section{summary}
%Now all summation goes in order of Fermi energy $E_F$, so we can replace closed-channel component $h_{2\vk}$ with $h_{2,0}=h_{2\,\vk=0}$
%\begin{equation}
%h_{1\vk}=\frac{\Delta}2\sqrt{\frac{(1-4h_{2\vk}^{2})}{(\epsilon^{ab}_{\vk}-2\mu+  h_{2\vk}^2\eta)^{2}+\Delta^{2}}}
%\end{equation}
%Gap equation 
%\begin{equation}\label{eq:20101004:gapS}
%\nth{{t_{0}}(\mu)}=\sum_{\vk}
%\br{\nth{\epsilon_{\vk}}-\frac{\sqrt{(1-4h_{2\vk}^{2})}}{\sqrt{{(\epsilon^{ab}_{\vk}-2\mu+  h_{2\vk}^2\eta)^{2}+\Delta^{2}}}}}
%\end{equation} 
%\begin{gather}
%{t_{0}}(\mu)=\br{1-\tilde{U}\tilde{ K}}^{-1}\tilde{U}\label{eq:20101004:t011}\\
%\tilde{U}_{\vk\vk'}=\nth{2} \br{U_{\vk\vk'}-\frac{Y_{\vk\vk''}Y_{\vk''\vk'}}{2(E^{0}-\eta+2\mu)}}\label{eq:20101004:tu11}\\
%{K}=\nth{\epsilon_{\vk}}\delta_{\vk\vk'}\label{eq:20101004:tk11}
%\end{gather}
%Here Eqs. (\ref{eq:20101004:t011}-\ref{eq:20101004:tk11}) follows the same two-body formula for zero-energy T-matrix element.  However, the detuning is shifted by a many-body quantity $\mu$ that should be determined by solving gap equation with number equation.  
%And number equation
%\begin{equation}\tag{\ref{eq:20100909:number}}
%N=\sum_{\vk} \left(1-\sgn_{k}\sqrt{1-4 h_{1\vk}^2-4 h_{2\vk}^2}\right)
%\end{equation} 
%
%\subsection{$G_0$ and $\eta$}
%%In the full region, $\eta$ should be close to the binding-energy of the closed-channel bound-state all the time.  The sweep in detuning should always be smaller than $\eta$ in order.  So maybe it is fine to take it as constant of binding energy $E^0$.  However, this quantiy is a two-bdoy quantity and very specific for different case. 
%%
%%$G_0=\alpha\phi_0$, where $\phi_0$ is normlized wave function of two-body closed-channel bound-state, strictly two-body.  Also, it is different from case to case. Its Pauli exclusion effect only shows up in many-body physics. It is determined by the real size of closed-channel bound-state and therefore relates to binding energy $E^0$.  
%%
%%For the same $a_s$ or $t_0$, $G_0$ or $\eta$ can take different value for different resonance.  Maybe we can take it as a different exogenous parameter. But we still need to find a way to estimate/deduce it from measurable quantities.  It is possible that they only relate to many-body quantities (same number for different densities maybe?)  because they are pertinent in many-body physics.
%%
%%\subsection{more on $G_{k}$ and $\eta$}
%$\phi^{0}$ the w.f. of two-body closed-channel bound state varies in different Feshbach resonance and is a case-specific quantity.  Its relevant property is included in the general two-body Feshbach resonance formula's $\Delta{B}$.  However, its other properties does matter in narrow resonance problem.  Fortunately, this can be simplified if it is a close-to-threshold bound-state, that mostly is outside the potential.  Therefore, most weight of $\phi^{0}$ follows the simple form 
%\[
%\phi^{0}(r)=C\frac{e^{-\kappa{r}}}{r}
%\]
% where $\frac{\hbar^{2}\kappa^{2}}{2m}=E^{0}$ and $\kappa>0$.  $C^{2}=\frac{\kappa}{2\pi}$.  Convert it into k-space
% \begin{equation}\label{eq:20101011:phi}
%\phi^{0}_k=\nth{(2\pi)^{\frac32}}\int{d^{3}\vr\phi^{0}(r)e^{-i\vk\cdot\vr}}=\nth{\pi}\frac{\kappa}{k^{2}+\kappa^{2}}
%\end{equation}
%Now we can estimated $G_{k}^{2}\eta$. $G_{k}=\alpha\phi^{0}_{k}$ where 
%\[|\alpha|^{2}<n\sim{E_{F}^{3/2}}\]
%\[
%\abs{\phi^{0}_{k=0}}^{2}=\nth{\pi^{2}\kappa^{3}}\sim\nth{(E^{0})^{3/2}}
%\]
%And $\eta\sim{E^{0}}$, So 
%\[
%G_{k}^{2}\eta<G_{k=0}^{2}\eta<E_{F}(\frac{E_{F}}{E^{0}})^{1/2}
%\]
%And we have $E_{F}\ll{E^{0}}$, so this term is much smaller than $\mu\sim{E_{F}}$ in BCS side. For now, we just omit this term.  Considering $G_{k}^{2}\ll1$, we can expand $\sqrt{1-4G_{k}^{2}}$ in \eef{eq:20101011:number} and \eef{eq:20101004:gapS} 
%\begin{gather}
%\nth{{t_{0}}(\mu)}=\sum_{\vk}
%\br{\nth{\epsilon_{\vk}}-\frac{1}{\sqrt{{(\epsilon_{\vk}-2\mu)^{2}+\Delta^{2}}}}
%+\frac{2\alpha^{2}|\phi^{0}_{k}|^{2}}{\sqrt{{(\epsilon_{\vk}-2\mu)^{2}+\Delta^{2}}}}}\label{eq:20101011:gapa}
%\\
%N=\sum{1-\frac{(\epsilon_{\vk}-2\mu)}{{\sqrt{{(\epsilon_{\vk}-2\mu)^{2}+\Delta^{2}}}}}
%+2\alpha^{2}|\phi^{0}_{k}|^{2}\frac{(\epsilon_{\vk}-2\mu)}{{\sqrt{{(\epsilon_{\vk}-2\mu)^{2}+\Delta^{2}}}}}}
%\label{eq:20101011:gapb}
%\end{gather}
%%Here $\phi^{0}_{k}$ takes the form in \eef{eq:20101011:phi} and we have an extra exogenous parameter $\kappa$ ($\approx{}E^{0}$) (which seems to independent of the normal Feshbach resonance formula).
%%\begin{equation}
%%a_{s}=a_{bg}(1+\frac{\Delta{B}}{B-B_{0}})
%%\end{equation}
%%This seems to make sense that in Feshbach resonance, what is really matter is the coupling energy, which involves  the integration of $Y_{kk'}$ and $\alpha^{2}\phi^{0}$ ($\alpha^{2}$ is the closed-channel weight in two-body physics).  The resonance is strong even when $\alpha^{2}$ is small but if the coupling strength $Y_{kk'}$ is strong.  However, the Pauli exclusion only cares about the wave function or $\alpha$, not coupling strength.   Therefore, we do need one extra parameter to describe this.  
%%
%%A more convenient choice for $\phi^{0}$ is to normalized to total number N, i.e. $\sum|\phi^{0}_{k}|^{2}=N$.  So $0<|\alpha|^{2}<1$ is the ratio of atoms in closed-channel to total number.  For this normalization, we find 
%%\begin{equation}
%%\phi^{0}_{k}=\sqrt{8\pi\kappa\rho}\nth{k^{2}+\kappa^{2}}
%%\end{equation}
%%
%% \cite{shizhongSumRule} derived  (eq. (25))
%%\begin{equation}\label{eq:20101025:fr}
%%F(\vr)=\frac{m\Delta}{4\pi\hbar^{2}}\frac{1-r/a_{s}}{r}
%%\end{equation}
%%
%%we have the pre-factor $\frac{m\Delta}{4\pi\hbar^{2}}$ for the normalization and with \cite{Leggett} (Eq. (4.A.17)), we can express the relation between $\Delta$ and $\alpha$.  (assume $\kappa$ varies slowly with detuning)
